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We have calculated the non-equilibrium quasiparticle and phonon distributions f(E), n(fl), where 
E and Q are the quasiparticle and phonon energies respectively, generated by the photons of the 
probe signal of a low temperature superconducting resonator SR operating well-below its transition 
temperature T c as the absorbed probe power per unit volume P a bs was changed. The calculations 
give insight into a rate equation estimate which suggests that the quasiparticle distributions can be 
driven far from the thermal equilibrium value for typical readout powers. From f(E) the driven 
quasiparticle number density N qp and lifetime i> were calculated. Using N qp we defined an effective 
temperature to describe the driven f(E). The lifetime was compared to the distribution averaged 
thermal lifetime at T^j and good agreement was found typically within a few percent. We used f(E) 
to model a representative SR. The complex conductivity and hence the frequency dependence of 
the experimentally measured forward scattering parameter S21 of the SR as a function of P a bs 
were found. The non-equilibrium S21 cannot be accurately modeled by a thermal distribution at 
an elevated temperature T£i having a higher quality-factor in all cases studied and for low P a bs 
T'21 ~ T^f. Using r r and N qp we determined the achievable Noise Equivalent Power of the resonator 
used as a detector as a function of P a bs- Simpler expressions for Tjv as a function of P a bs were 
derived which give a very good account of and also N qp and r r . We conclude that multiple 
photon absorption from the probe increases the quasiparticle number above the thermal background 
and ultimately limits the achievable NEP of the resonator. 
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I. INTRODUCTION 



Superconducting resonators (SRs) with high quality 
factor Q operating at low reduced temperatures T/T c ~ 
0.1, where T c is the superconducting transition tem- 
perature, are used not only as ultra-sensitive detec- 
tors of individual quanta or incident power for applica- 
tions in sub-millimeter, millimeter, optical, X- and 7-ray 
astrophysics fi _ — but also as elements of Qubits for quan- 
tum computing, They are also needed as elements of 
microwave Superconducting Quantum Interference De- 
vice (SQUID) multiplexers^ and also some SQUID ge- 
ometries more directly. Despite this technological im- 
portance, we are unaware of any detailed analysis of 
the effect of the interaction of a flux of microwave pho- 
tons of frequency v p <C 2A(T)/h, where 2A(T) is the 
temperature-dependent superconducting energy gap and 
h is Planck's constant, on the superconducting state at 
low temperatures. One might assume that if hv p < 2A 
the photon interaction cannot change the quasiparticle 
number since the photon cannot break a Cooper pair. 
Whilst true at the single quantum level, that assump- 
tion ignores the effect of a flux of photons comprising 
very many quanta which might be used in a typical 
experiment and the competing effects of multiple pho- 
ton absorption and quasiparticle scattering and recom- 
bination on the driven quasiparticle distribution f(E) 
where E is the energy. Understanding the effect of non- 



equilibrium quasiparticles on a Qubit is certainly a topic 
of current interest; 10 ! 11 where non-equilibrium quasipar- 
ticles may be a limiting factor on Qubit energy relaxation 
times. Our interest is particularly in the context of non- 
equilibria in SRs used as quantum and power detectors 
although we would emphasize that our calculations apply 
to SRs in general and the mechanisms and solutions we 
describe are common to all SR applications. Indeed the 
results presented below are independent of the particular 
geometry or application provided the effect of geometry 
on the power absorbtion is considered. 

In the context of power detection, SRs are sensitive 
to changes in detected power because of the dependence 
of the surface impedance Z s of the superconductor on 
quasiparticle density and the dependence follows from 
the complex conductivity a described by Mattis and 
BardeenjH The SR is embedded in an electrical read- 
out circuit and is driven by a microwave probe signal of 
frequency v v <C 2A(T)/h close to the circuit resonant 
frequency i>q . If the quasiparticle density is changed, for 
example by absorption of a high frequency photon sig- 
nal of sufficient energy to break Cooper pairs hv$, > 2A, 
where v<& is the detected photon frequency, Z s is changed 
and the change can be monitored by measuring the 
change in the (complex) resonance transmission charac- 
teristic S^iiy) of the probe signal. The change relaxes 
back to the unperturbed state as energy is exchanged 
between the quasiparticles and the phonons of the su- 
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perconductor and ultimately the substrate. In this way 
very sensitive power detectors can be made. For exam- 
ple Ref. (QJ) estimated a dark Noise Equivalent Power 
(NEP), i.e. ignoring the achieved signal detection ef- 
ficiency, of 2 x 10~ 19 W/VHz at lowest readout power 
which was accounted for in terms of the generation- 
recombination noise of the quasiparticles and a limiting 
lifetime of 3 — 4 ms. It has been predicted that NEP's ap- 
proaching 10 -20 W/\/Hz or lower may be achieved using 
SRs. The readout naturally lends itself to frequency- 
division multiplexing where a large number of SRs each 
operating with a slightly different vq are coupled to a 
through-transmission line. A high Q resonator can be 
formed by lithographically patterning a low-T c supercon- 
ducting thin film such as Al (T c = 1.2 K) on a dielectric 
substrate. The substrate is held at the bath temperature 
T = T&. A number of SR geometries are possible in- 
cluding ring, half- and quarter-wave or lumped-element 
designs. 

In a recent paper we discussed the effect of the probe 
power on the resonator characteristic?^ The key point 
is that the readout is dissipative^ i 13 ' 15 and the absorbed 
power can be calculated knowing the embedding electri- 
cal circuit. Probe photons are absorbed by the quasi- 
particles of the SR which changes f(E), although the 
resulting distribution has not yet been calculated. In a 
real device other mechanisms may contribute to the dis- 
sipation for example dielectric or radiative losses! 2 - Here 
we focus on the dissipation associated with the real part 
of a. The aim of this work is to calculate the effect of 
a microwave drive at low temperatures (Tj/T c = 0.1) on 
the static, non-equilibrium quasiparticle and phonon dis- 
tributions in a superconductor as the probe-power levels 
are changed. We also derive simpler analytical expres- 
sions which give a good approximation to the key results 
of the full calculation. 

Energy relaxation processes of quasiparticles in a su- 
perconductor that couple to phonons comprise scattering 
of quasiparticles with absorption or emission of phonons, 
and pair-breaking and recombination of quasiparticles 
with absorption or emission of phonons of energy > 
2A. Energy escapes from the superconductor as phonons 
enter the substrate. The coupled kinetic equations that 
describe these interacting subsystems were derived by 
Bardeen, Rickayzen and Tewordlji^ and discussed in de- 
tail by Chang and Scalapino i 17 ' 18 In Ref. [lj] the coupled 
equations were linearized and solved for a variety of drive 
sources including microwaves. In Ref. [TH full non-linear 
solutions were obtained. Crucially however those solu- 
tions were obtained close to T c where microwave drive can 
lead to gap-enhancement effects. The kinetic equations 
have been used to investigate the non-linear effect of high 
energy photon interactions at low temperatures*^— To 
our knowledge no full non-linear calculations exist of the 
effect of microwave drive at low temperatures on SRs. 

The structure of this paper is as follows. In Section HTl 
we give an estimate of the power densities where non- 
equilibrium effects are likely to become important in a 



SR. In Section IIIII we discuss the general properties of 
the coupled kinetic equations and derive the form for the 
drive and subsystem power-flow terms necessary to en- 
sure energy conservation whilst also discussing the effects 
of out-diffusion in a real geometry and the method used 
to calculate distribution averaged recombination times. 
In Sect ion HVl we give numerical parameters used in calcu- 
lations to describe a clean thin- film superconductor (Al), 
and in IIV Al we describe the model used to calculate S21 
for a representative SR. Section [V] describes the numer- 
ical method. In Section I VII we show solutions for the 
non-equilibrium quasiparticle and phonon distributions 
of a driven SR operating with Tj,/T c = 0.1 as a function 
of P a bs including calculations of driven quasiparticle den- 
sity N qp , an effective temperature determined from 
N qp , and of the distribution- averaged relaxation time T r . 
We calculate the driven for a representative Al SR 
under the same powers. We also estimate the effect on 
the achievable NEP for a quantum SR detector using 
these results. Section IVTIl describes two models both giv- 
ing a reasonable but simpler account of the results. In 
Section lVIIII we discuss the implications of the work with 
concluding remarks. 

II. EQUILIBRIUM ESTIMATE 

In this section we use the equilibrium interaction times 
derived by Kaplan et al££ to estimate the power den- 
sities where non-equilibrium effects are likely to occur 
in a SR at low temperatures. Ref. gives expres- 
sions for the thermal equilibrium lifetimes: t s {E,T) for 
the scattering of quasiparticles of energy E, r r (E,T) 
for recombination to Cooper pairs, T0 S (O, T) for the 
scattering of phonons of energy f2 and T P &(Q, T) for a 
phonon to break a pair. These are calculated in terms 
of characteristic quasiparticle and phonon lifetimes To 
and Tq . The volume density of thermal quasiparti- 
cles is N(T) = 4AT(0) /~ p(E,A(T))f(E,T)dE, where 

p(E,A) = E/\/E 2 — A 2 is the normalized quasiparticle 
density of states, f(E, T) = (exp(E/k b T) + 1) _1 is the 
thermal Fermi distribution at temperature T (in contrast 
to f(E) the driven distribution), fc& is Boltzmann's con- 
stant and N(0) is the single-spin density of states at the 
Fermi energy. Consider the case of an Al resonator with 
a probe signal of v v ~ 4 GHz. Since hv v ~ 2A/20 it is 
not clear that the probe-signal is capable of changing the 
equilibrium N(T) and hence A(T). Our measurements 
and modeling show that typical experimental read-out 
powers dissipate of order 1 — ^ 100 aW//jm 3 in an Al SR 
used as a power detector .^±1 

The probe photons are absorbed by the quasiparti- 
cles changing their energy distribution f{E). We ex- 
pect the appearance of at least one peak in f(E) around 
E = hv p due to absorption of the monochromatic probe 
photons by the large density of quasiparticles near to the 
gap, perhaps more peaks if P a i> s is sufficiently high. As- 
sume that energy relaxation of an excited quasiparticle 
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FIG. 1: (Color online) Quasiparticle energy E$ as a function 
of P a bs for an Al film at 3 bath temperatures assuming cooling 
by quasiparticle-phonon scattering alone. 



can only occur by scattering with emission of phonons 
and all of the emitted phonons are lost from the film 
with an energy-independent time r;. Indeed, for thermal 
distributions at low reduced temperatures quasiparticle- 
scattering times are significantly shorter than recombi- 
nation times' However if f(E) has any non-equilibrium 
quasiparticles above E = 3A the phonon emitted in 
scattering may have SI > 2 A which can break a pair. 
This is the onset of non-equilibrium effects since the 
number of quasiparticles can be changed. We can ob- 
tain a naive estimate of the probe-power level for this 
to occur by assuming that all of the available quasi- 
particles at T b i.e. N(T) = N(T b ), are driven to the 



same unknown energy Eg and that scattering is the 
only energy loss mechanism for the quasiparticles, so 
that P abs = E 5 N(T b )/T s {E s ,T b ). Figure CD shows E 5 
for three bath temperatures as a function of P abs calcu- 
lated for an Al film with A(0) = 180 fieV, T c = 1.17 K, 
AT(0) = 1.74 x 10 4 ^eV" 1 ^nT 3 , and t = 438 ns. At a 
typical operating temperature T b /T c = 0.1 all the quasi- 
particles may have Eg > 3 A. Even at T b /T c = 0.2 E$ is 
close to the pair-breaking threshold due to reabsorption 
of the scattered phonon for typical SR drive powers. We 
emphasize that this estimate is very conservative because 
non-equilibrium effects occur as soon as there are any 
non-thermal quasiparticles with E > 3A and the degree 
of non-equilibrium depends on the probability of pair- 
breaking before loss by the scattered phonons, Ti/T pb . 
At low temperatures r pb — Tq. For a 100 nm Al film 
on sapphire or Si we estimate that t; ~ r pb 2^ This sug- 
gests that the probe-signal may break Cooper pairs in the 
driven SR even for very low P abs and scattering without 
pair-breaking may not be a sufficient energy-relaxation 
mechanism. Then a full non-equilibrium description is 
required and we will not be surprised if the power densi- 
ties to observe non-equilibrium effects may be lower than 
suggested by Fig. [T] 



III. NON-EQUILIBRIUM RESONATORS 

In this and subsequent Sections we describe and use 
a full non-linear solution of the kinetic equations for the 
coupled quasiparticle and phonon systems. The rates of 
change of the distribution functions, f(E) for the quasi- 
particles and n(Sl) for the phonons, are given in Eqs. (7) 
and (8) of Ref. 18. Substituting tq and Tq these become 



ffi^ = i a JE, u„) — ; 1 r „ f dsisi 2 P (E + si) ( 1 — - 

dt qpK ' pJ r {k b T c f J n ; V E{E + VL) 

(/ (E) [1 -f(E + n)] n (fi) - [1 - / (E)] f (E + (I) [n (O) + 1]) 

1 pE—A , A 2 



■ J dnn 2 P (E - fi) ^i - 



MhT c ) 3 J KV J \ E(E-tt)J (1) 

(/ (E) [l-f(E- Q)] [n (fi) + 1] - [1 - / (E)] f (E - SI) n (Si)) 
1 f°° „ _ „ / A 2 



/ dSlSl 2 p(Sl - E) [ 1 
Je+a V 



T {k b T c f J E+A rv '\ E(Sl-E) / 
(/ (E) f(Sl-E) [n (SI) + !]-[!-/ (E)} [I - f (SI - E)} n (Si)) , 
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and 



J dEp{E)p{E + SI) ^1 



dn(Sl) 

dt ttt$A(0)Ja — rv-"-v- ■ £(£7 + fi)_ 

(/ (£) [1 - / (E + Si)} n (Si) - [1 - / (£)] / (£7 + SI) [n (SI) + 1]) 
fi-A / 



dEp(E)p{Sl - E) 1 + — — — x (2) 



ttt*A(0) J a V £(ft-£7), 

([1 - / (£7)] [1 - / (fi - £0] n (fi) - / (£7) / (n - £7) [n (Q) + 1]) 
n(fi) - n(fi, T fe ) 

) 



where n(f2,Tb) = (exp(f2/fcf,Ti>) — 1) is the Bose dis- equilibrium f(E). In the static, driven situation to be 
tribution evaluated at T b . The term I qp (E,v p ) is the solved df(E)/dt = dn(Sl)/dt = 0. 
source term due to the photons at energy E and quan- 
tifies the drive of the microwave probe. The energy gap 
is modified from its equilibrium value and is determined 
self-consistently so that 

1 _ f 1 ~ 2/(£7) , . A. Microwave drive term 

N(0)V B cs ~ J a VE 2 - A 2 ' 1 ' 

where Vbcs is the BCS interaction parameter. Eq. ([3]) The form of I qp (E, hv v ) was calculated by Eliashberg 24 

calculates the non-equilibrium A using the non- and 



I qp {E,v p ) = 2B 



p(E + hv p , A) 
p(E-hv p ,A) 



A 2 



where the rate coefficient B needs to be determined for 
low temperatures. A third term arises in Eq. (jj) if 
hv p > 2A , which is not considered here since we are 
investigating the effect of sub-gap photons. Eq. H] de- 
scribes both absorption and emission of single photons. 
Ref. [It] used a different rate coefficient to describe the 
interaction of microwaves at normal incidence to a su- 
perconducting film close to T c which, in our notation, 
would be B' oc H 2 R n / dN {Q)v 2 with H the magnetic field 
strength, R n the normal-state square sheet resistance and 
d the film thickness . For a SR at low temperatures 
the field does not interact with a normal-state metal in 
that geometry, neither is the penetration of the field into 
the superconductor determined by the normal-state pa- 
rameter. We take a different approach which emerges 
naturally by considering energy conservation. Assuming 
uniform absorption, the power absorbed per unit volume 



E(E + hu p ) 
A 2 

E (E - hv p ) 
I 



[/ (E + hv p ) - f (E)} 
[/ (E) -f(E- hup)] 



(4) 



of the resonator is 



P abs =4N{0) dEI qp (E,v p )Ep(E,A). (5) 

J A 

We solve Eqs. [I] and [5] numerically so that writing 
I qp (E,i> p ) — BK qp (E^v p ) we can also include a power 
absorption error term 

5P = 4N(0)B / dEK qp (E, v p )Ep(E, A) - P abs , (6) 

J A 

where B needs to be determined to satisfy this equation, 
and the static non-equilibrium solution sought is SP = 0. 



B. Quasiparticle-phonon power 

The power flow from the quasiparticles to the phonons 
can be found by recognizing that the sum of the inte- 
grals on the right-hand-side of Eq. (fTJ) gives the total 
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rate of change of f(E) due to interactions with phonons, 
Iqp-c/>{E). The energy leaving the quasiparticles per unit 
volume per unit time is 

poo 

Pgp-^ANiO) dEI qp _^(E)Ep(E,A). (7) 

J A 

We define the fractional quasiparticle-phonon power flow 
error term 

C Pahs — Pqp~4> / q\ 

*abs 

C. Phonon cooling term 

Energy is lost from a SR as non-equilibrium phonons 
are lost into the substrate. The energy leaving the 
phonons per unit volume of the film per unit time is 
given by 

P ,_ b = r^Mn n(0) ~ n(0,Tb) . ( 9 ) 

* Jo n 

With a Debye model the density of states is given by 
D{uj) = 3il 2 /fl^, is the Debye energy, and the sum 
over the phonon branches introduces an additional factor 
of three. We define a further error term 

b_ b = Pabs p p f- b . (10) 

*abs 

Eqs. and ()10|) provide an important monitor of the 
accuracy of the numerical solutions to the coupled equa- 
tions. 



D. Power absorption and the effect of geometry 

Power absorption in a resonator is dependent on its 
geometry. In modeling, we assumed that probe photons 
were uniformly absorbed in the SR and we ignored pos- 
sible out-diffusion of excited quasiparticles. We assumed 
Pabs was known. In practice geometric effects need to 
accounted for but these depend on the particular real- 
ization. For example a A/4-resonator is a useful geom- 
etry into which to couple an external pair-breaking sig- 
nal. The current density distribution in this device is 
described by J(x) — J(xo) sui(t:x/2xo), which is a max- 
imum J(xq) at the shorted end x = xq, and the peak 
power density is a factor close to 2-times higher than the 
average. For the same reason, as a detector the A/4-SR 
is most sensitive to changes in the quasiparticle density 
at its shorted end. Mirror currents in the ground-plane 
mean that the effective volume for the power absorption 
is upto twice that of the central conductor. Moreover in 
a A/4 resonator used as a detector out-diffusion of the 
excess quasiparticles generated by the in-coming signal 
must be minimized in order to maximize the detection 
sensitivity, which may be achieved for example by using 



a higher energy gap contact to the electrical ground. In 
practice other resonator geometries are used and the ef- 
fect of geometry and out-diffusion on power absorption 
can in principle be calculated. 

E. Recombination times 

In what follows we calculated distribution-averaged 
quasiparticle recombination times r r for the driven f(E). 
We used Ref. [13 Eq. (A9) to find the rate coefficient R 
averaged over f(E) and set r r = (2RN qp )~ 1 with N qp = 
4iV(0) f. dEp(E)f(E) the non-equilibrium quasiparti- 
cle density. If the detected power is small compared to 
the probe power this is the appropriate measure of the 
small-signal relaxation time. For a thermal distribution 
we find T r = (r r (T))g P , the distribution-averaged recom- 
bination time described by Kaplan et ai. 

IV. NUMERICAL PARAMETERS 

To describe the resonator we used material parame- 
ters appropriate for Al with N(0)Vbcs — 0.167 giving 
A(0) = 180 peV, T c = 1.17 K and we set T/T b = 0.1. 
We used N(0) = 1.74 x 10 4 peV~ 1 pm~ 3 , r = 438 ns and 
Tq' 1 = 0.26 ns22. The latter was calculated assuming that 
the appropriate value for the phonon density of states 
in the calculation of Tq h is a 2 D ap pro priate for a Debye 
model as given in Table II of Ref. |22|, an approach sug- 
gested in Ref. fl9L To be precise, in our view the param- 
eters needed as inputs for modeling are not collectively 
and with sufficient precision known from measurement or 
theory. Considering the pre-factors in Eqs. (TTJ), 0, 
and ([9]) , we found that the numerical inputs must satisfy 

2nN{Q)4 A tf D = i 

QN ion T (k b Tcf 

to allow a self-consistent solution where the power errors 
Eqs. dSJ) and (fit))) converged to zero. 

A. Parameters for a A/4-resonator 

In calculations discussed later and shown in Fig. [7] we 
investigated the effect of the driven f(E) on a represen- 
tative device modeling a A/4 microstrip resonator as in 
Ref. IbH . The complex conductivity a which is propor- 
tional to the normal-state conductivity <tjv was calcu- 
lated from Eqs. (3.9) and (3.10) of Ref. Q but using 
the non-equilibrium f(E). The surface impedance was 
calculated from a hence the propagation constant and 
the characteristic impedance of the SR. The modeled SR 
had a length of 7.6 mm, width 3 pm, film thickness of 
200 nm, dielectric thickness of 200 nm with e r = 3.8 and 
a saturation quality factor of 10 7 . We set the coupling ca- 
pacitance to be 5 fF and (Jn = 1-25 x 10 8 (Q m) _1 , which 



6 



would be typical for a clean Al film at low temperatures. 
We calculated a resonant frequency Uq = 3.934 21 GHz 
(hvo = 16.2 /ieV) with zero absorbed power. 



V. NUMERICAL METHOD 

A non-equilibrium solution for f(E) and n(Q) requires 
simultaneous solutions of Eqs. (QJ, ©, ©, and ([5]). The 
task is complicated since, for example, Eq. ([T|) contains 
terms such as f(E± Q) and likewise © requires knowl- 
edge of f(E± hi/p), and there are likely to be peaks in 
the driven distributions arising from the high density of 
states of quasiparticle near E = A. 

Eqs. ([J), ([2]), and ([5]) were solved using Newton's 
method. We discretized the distributions f(E) and n(fi) 
using a 1 /ieV grid with Ei = A + i — 1, = i 
and i 6 1 . . . N with N = 1000 so that that quasi- 
particle states up to ~ 6.5A(0) are considered. We 
formed the state vector a = [fi,B,n,i] T where T de- 
notes the transpose. We formed the error vector e = 
[dfi/dt,SP,drii/di\ T . The iterative procedure seeks to 
find e l+1 — using ot l+1 = a 1 — \ e ( a where 

the matrix J = dej/dak is the Jacobian of the partial 
derivatives and j, k G 1 . . . 2N + 2. It is possible to de- 
rive analytical expression for the derivatives making the 
Jacobian efficient to evaluate. The superscript I denotes 
the iteration number, x < 1 is a convergence parame- 
ter and we find 0.8 < x < 0.95 gives reasonably rapid 
convergence typically within 10 iterations. 

We assumed a starting thermal f° = f(Ei, T start ) with 
an initial temperature T star t ~ 2Tb The value chosen 
for T s tart did not affect the solutions obtained merely 
the number of iterations required to converge sufficiently. 
Using an earlier estimate of the non-equilibrium distri- 
butions reduces the number of iterations (or increases 
the precision for the same computation time) and is a 
useful approach if parameters such as power or phonon- 
trapping factors are being varied systematically. We 
chose n° — n(fij,Tj,) so that the phonons are initially 
at the bath temperature. The aim was to find f l (E) and 
n'(fi) such that |£g P _4,|,l4_bl < 0.1% for both power 
transfer error terms. All solutions shown below exceed 
these convergence criteria in some cases by nearly an or- 
der of magnitude. 

We found that in our solutions the non-equilibrium gap 
calculated with Eq. ^ changed very little from A(T&) 
and by a maximum of i5A ~ 50 nV. For this reason 
we did not allow A to change in the simulations. For the 
microwave drive we restricted v v to match the discretized 
distributions. This means that the onset of any photon 
induced peaks occurs in well-defined bins of f(E{) and 
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FIG. 2: (Color online) Non-equlibrium distribution for ab- 
sorbed power 2 FW/nm 3 with T b /T c — 0.1 and ri/r pb — 1. 
The continuous curve is a Fermi distribution f(E, T^) hav- 
ing the same quasiparticle density. The inset shows semi- 
log plots for powers of (a) 2 FW/fim 3 (full), (b) 20 aW/^m 3 
(dashed) and (c) 0.2 aW//im 3 (dot-dash). The associated 
straight lines, with the same line styles, show the thermal 
distributions having the same number density of quasiparti- 
cles as the driven f(E) in each case. 



VI. SOLUTIONS FOR NON-EQUILIBRIUM 
RESONATORS 

In this section we show results of the modeling. Fig- 
ureHJshows f(E) with v p = 3.880 4 GHz (hv p = 16 /xeV), 
which is close to v a for the A/4-A1 SR that we later use 
as an example, with T/T c — 0.1, P a b s = 2fW//xm 3 
and Ti/Tpb = 1. The multi-peaked structure is consis- 
tent with sequential single photon absorption; the drive 
term of Eq. (|4]) only describes single photon events. The 
occurrence of this multi-peaked structure is expected if 
the SR is driven far from equilibrium. Pleasingly this 
structure emerges in the very first iteration of the nu- 
merical method. In addition physically unrealistic dis- 
tributions (where for example f(E— hv) < f(E) at 
E < 3A) were never found. Figure [5] includes a ther- 
mal distribution f(E,T^) where is defined so that 
47V(0) jf^° f{E,T* N )p{E)dE = N qp i.e. the thermal dis- 
tribution having the same number density of quasiparti- 
cles as the driven case. The non-equilibrium nature of 
f(E) becomes further apparent in the semi- log plots in 
the inset which are calculated for P aos of (a) 2 fW /fim 3 , 
(b) 20 aW//im 3 and (c) 0.2 aW//im 3 where the presence 
of quasiparticles with E > 3A also showing multiple pho- 
ton induced structure can be seen. These quasiparticles 
arise from absorption of 2A-phonons by quasiparticles. 
A further much-reduced feature (not plotted) at E > 5A 
is also found. A recurring feature of these solutions is 
that the driven low-energy f(E) shows excess densities 
of quasiparticles at energies of order E < A+10hv p above 
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FIG. 3: (Color online) Contributions to the phonon power 
flow for absorbed power 2 fW / /im 3 . The inset shows semi-log 
plots for absorbed powers of (a) 2 fW//im 3 , (b) 20 aW//im 3 
and (c) 0.2 aW//im 3 . all with Ti/T p b = 1 



FIG. 4: (Color online) Contributions to the phonon power 
to the bath P(fl)^,- b and inset the associated f(E) for P abs = 
50 aW/^m 3 , T b /T c = 0.1, for two values of ri/r pb = 0.5 (full 
lines) and 2 (dash-dot lines) with Ti/r pb = 1. 



the equivalent distributions (the dashed lines). The 
effect of this distortion is to increase the power carried by 
low energy phonons Q < 2A by scattering compared to a 
thermal distribution and these phonons are more easily 
lost from the SR providing an efficient cooling mecha- 
nism. The distortion from the equivalent thermal distri- 
bution increases as P a & s increases, as does the number 
of photon peaks. At energies E ~ 3A — hhv v the cal- 
culated f(E) is increased from the by-eye straight line 
and at the same time for E > 3A it is indeed the case 
that f(E — hv) < f(E) for some E. The magnitude of 
both effects are power dependent, which arises from the 
competing contributions of non-equilibrium 2A-phonon 
re-absorption and the stimulated emission of photons in- 
herent in Eq. (g]). 

Figure [3] shows the corresponding contributions to the 
phonon-bath power flow integral of Eq. ((9]), P(fi)<£_(„ 
for the same drive conditions where the presence of non- 
equilibrium 2A-phonons is seen. An additional feature at 
n = 4A is also found. There are two distinct contribu- 
tions to P(fl)^-b- At low energies we see phonons arising 
from the scattering of low energy quasiparticles towards 
the gap, and there is structure consistent with the peaks 
in f{E). Structure on the low energy side of the phonon 
peaks is also seen which is expected as the driven f(E) 
scatters to lower energies and the rate of this scattering is 
reduced by the occupation of the final states by the driven 
distribution itself, despite the increasing phonon density 
of states that would be available for the scattering to oc- 
cur. At SI > 2A we see a second distinct contribution 
to P(tt)<p- b . This power is transferred to the substrate 
by pair-breaking phonons with SI > 2A, which are gener- 
ated not just by recombination of the excess f(E) itself 
but also by the scattering (and recombination) of those 
quasiparticles with E > 3A. To quantify the fraction 
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FIG. 5: (Color online) The effect of P abs = 50 aW/pm 3 as 
a function of Ti/r pb : (a) Quasiparticle density, (b) Effective 
temperature T^r, (c) Fraction of power carried by 2A-phonons 
and (d) Recombination time for the non-equilibrium f(E). 



of the power carried by phonons with SI > 2A we define 
V2A = f£d£l(n(Sl)-n(Sl,T b ))/ J °° dfi(n(fi) -n(0, T b )). 
In the main plot of Fig. [31 7/2A = 0.16. Figure @] shows 
contributions to P(Q)<f)-b for two values of Tijr v b- Some- 
what counter-intuitively increasing r;/r p b increases the 
contribution at the lowest SI to P(0) 9 j_j, whilst simulta- 
neously increasing the contribution from pair-breaking 
phonons SI > 2A. The effect on f(E) is evident in 
the inset. As n/T p b is increased more 3 A quasiparticles 
are generated, these in turn generate more pair-breaking 
phonons which can be reabsorbed before being lost from 
the film. Figure [S] shows the effect of changing Ti/r p b 
772 A, and T r for P a b s = 50 aW/^m 3 . In- 



on N qp , Tfr, 
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FIG. 6: (Color online) The effect of absorbed power with 
ti/t p i, = 1: (a) Quasiparticle density, (b) Effective temper- 
ature Tff, (c) Fraction of power carried by 2A-phonons and 
(d) Recombination time for the non-equilibrium f(E). The 
full curve in (b) is an anayltical expression. The full curve in 
(d) is the distribution-averaged thermal recombination time 
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FIG. 7: (Color online) The effect of absorbed power on S21 
(full lines) with Ti/rpt, = 1 and (dashed lines) the transmis- 
sion calculated for the nearest thermal distribution f(E, TJi) 
giving the same resonant frequency: (a) P a t, s = at Tt, 
(b) Pabs = 0.1 aW//im 3 , (c) P abs = 2 aW//im 3 and (d) 
Pabs = 50 fW //im 3 . 8v is referenced to vq with P a b 3 = 0. 



creasing Ti/r p b increases both N qp and while r r is 
reduced. In combination this increases the relative con- 
tribution of the recombination phonons to the power 
flow to the bath despite the increasing probability of re- 
absorption, and for this reason 772 a increases somewhat. 
We find that for fixed power the solutions are related by 
N q p/T r (l + Ti/Tpi,) = k where k is a constant independent 
of 772A although we emphasize that both N qp and r r are 
driven non-equilibrium values. 

Figure [6] shows the effect at T/T c = 0.1 of varying 
P a bs on N qp , T^, 772 a and r r . Figure [Bfa) shows that, 
for all P a bs, N qp exceeds the undriven thermal density 
calculated at T and, to emphasize, for Tb/T c = 0.1 we 
calculate N qp = 0.1 fim~ 3 . Readout power significantly 
changes the driven, static N qp . Figure EJb) shows that 
is enhanced above Tj, for all read-out powers studied. 
The full curve shown is an analytical expression described 
later in Sec. I VIII The distortion of f(E) from even the 
nearest thermal distribution as P a bs is increased means 
that 772A shown in Fig.JSJc) is also a function of P a bs- At 
the lowest powers studied P a bs ~ 0.1 aW//im~ 3 much of 
the power leaving the film is carried by recombination 
phonons, which is as expected given our earlier estimate 
showing the inefficiency of scattering in the energy re- 
laxation. As P a bs is increased ?72A is reduced and more 
power is carried by < 2A-phonons emitted by quasi- 
particle scattering. For P a bs — > we find 772A — > 0.6 
which interestingly is the result found for high-energy 
interactions hv 3> 2A. Figure EJd) shows that r r is re- 
duced as P a b s is increased and, as expected, mirrors the 
increase in N qp . The full curve in (d) is the distribution 
averaged thermal recombination time at T^, (T r (T^)) qp , 



which gives a very good description of the recombination 
time of the driven f(E) typically within a few %. The de- 
viation increases with P a bs which again is expected: the 
distortion of f(E) with P a b s from the distribution in- 
creases the available final phonon densities of states for 
recombination. 

Figure [7] shows (full lines) calculated IS21I 2 f° r the 
driven f(E) for absorbed powers: (a) P a bs — at Tb, 
(b) Pabs = 1 aW//jm 3 , (c) P a bs = 2 fW/Vm 3 and (d) 
Pabs = 50 aW/fim 3 . The dotted lines show l^il 2 for the 
nearest thermal distribution /(E^T^) giving the same 
vq in each case. A general characteristic of all calculated 
transmission curves is that the driven S21 has a higher 
Q (it is deeper and narrower) than the nearest T 2 *i P re ~ 
diction. This arises due to occupation of final states for 
absorption by the driven distributions. For increasing 
Pabs the divergence increases. This is a further effect of 
the increasing distortion of f(E) as a function of P a bs 
observed in relation to Fig. [2] The inset shows that 
gives a reasonable account of for the range otP a bs 
considered, particularly at low powers. 

The intrinsic limiting NEP of a superconducting de- 
tector in thermal equilibrium is determined by the ran- 
dom generation and recombination of quasiparticles^ In 
thermal equilibrium NEPq-r = 2 A/77 \J N qp V/ t* where 
77 is the fraction of detected power Pdet coupling to the 
quasiparticles, V is the volume of the SR, and the ef- 
fective recombination time r* = r r [1 + t//t po ] /2. Here 
the factor of two arises because, as noted by Kaplan 
et at and others, r r is that of a single quasiparticle 
whereas two quasiparticles are lost in each recombina- 
tion event i 22 ' 26 ' 27 For detection of high energy photons 
hi/g, 2A, 77 ~ 0.59*i2~— . If we assume that Pdet is small 
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FIG. 8: (Color online) The effect of absorbed power on 
generation-recombination limited NEP for 2 values of signal 



detection efficiency w, resonator volume V 
ri/r pb = 1. 



1000 /jm J and 



compared to P a bs so that 5N qp /N qp is small (as it must be 
for a linear detector) then the relevant N qp and 
already calculated. Figure [3] shows NEPq-r as a func- 
tion of P a bs for a SR with V = 1000 /mi 3 and t;/t p & = 1 
for 2 values of rj, and we have assumed that the equilib- 
rium expression for NEPg-r applies for the driven case. 
Pabs determines the limiting NEP and even for the lowest 
P a bs studied, NEP G _ R - 1 - 5 x 10- 19 W/VHz. 



VII. ANALYTICAL POWER MODEL 



The preceding calculations provide insight into the ef- 
fects of P a bs with Up <C 2A/h at low reduced temper- 
atures on low-T c SRs. However in many situations, for 
example for estimates of performance or for extrapola- 
tion to other materials, an expression to approximate the 
key results would be extremely powerful. Recombina- 
tion determines the overall time-evolution of the driven 
system, even though we have shown that only a power- 
dependent fraction 772A of P a bs is carried by the recom- 
bination phonons. An approximate equation giving an 
estimate of as a function of P a bs can be derived con- 
sidering energy conservation so that 



^ dPrj2A= [ T dT ^ BCS ^ 

" J T b ( T r{T)*) qp 



(12) 



where the denominator on the right-hand side is the 
distribution-averaged effective thermal recombination 
time and t*(T) = r r {T) [1 + Ti/r pb ] /2. C BC s is the 
BCS specific heat capacity which comprises two terms' 
The first is the quasiparticle heat capacity C qp = 
AN(<S)d/dT (/~ dEEf(E)p(E)) and the second arises 
because the quasiparticle energies E themselves change 



due to their dependence on A. At the (effective) tem- 
peratures considered here dA/dT ~ 0. The data of 
Fig. [6jc) were fitted to a log-linear model giving 772A = 
-0.031n(P Qf , s /aW/im- 3 ) + 0.384. At very low absorbed 
powers P — > we find 772 a — > 0.6. With the same limit 
Eq. (|T2"j) was solved. The result is shown as the full line 
in Fig. |5] (b). The account of as a function of P a b s is 
very satisfactory. 

We found that the functional form of Tjy with P a bs can 
be further approximated by the simpler expression 



1 



V2A(P abs ) 

x T^exp 



1 



1 + n/Tpb 
-2A(T* N ) 



k h T, 



b 1 N 



T exp 



-2A(7V 
hTb 



(13) 



where f72A(p a(J3 ) is the fraction of power carried by 2A- 
phonons at P a bs ■ For the Al film modeled here we found 
S s = 3.4 x 10 10 W/m 3 K. This function is indistinguish- 
able from the full curve plotted in Fig. [6fb) and gives a 
good account of the effect of P a bs ° n for the param- 
eter space studied. Both Eqs. (TT21 and (fll?j) provide a 
straightforward route to estimate N qp and r r as a func- 
tion Of Pabs- 



VIII. DISCUSSION AND CONCLUSIONS 

We have calculated the non-equilibrium distributions 
of quasiparticles and phonons, f(E), n(f2) generated by 
a flux of low-energy photons hv p <C 2A as a function 
of P a bs for a thin-film superconducting resonator at low 
temperatures T/T c = 0.1. Driven f(E) deviate from 
thermal-like distributions exhibiting structures associ- 
ated with multiple probe-photon absorption and emis- 
sion for all P a b s studied. All calculated n(f2) show pair- 
breaking phonons ft > 2A for all P a b s studied. The den- 
sity of driven quasiparticles exceeds their thermal density 
at the bath temperature, confirming a simple estimate 
based on energy conservation using thermal scattering 
times. The driven f(E) can be characterized in terms of 
an effective temperature which also gives a good ac- 
count of the distribution averaged, driven recombination 
time r r and this can be very- well approximated using sim- 
pler expressions to calculate the thermal recombination 
time at T^. Using N qp and r r a (dark) detector Noise 
Equivalent Power can be calculated. We find that dissi- 
pation limits the achievable NEP in the range of P a bs con- 
sidered indicating a minimum NEP ~ lx 10~ 19 W/VHz 
although we emphasize this depends on the actual ab- 
sorbed power. 

Ref. [l3| measured a A/2 Al resonator where we ex- 
pect out-diffusion should be minimized. That work es- 
timated a limiting effective quasiparticle temperature of 
order 160 mK with T b = 100 mK. N qp - 20 - 70 /Jm" 3 , 
r r 3.5 — 0.5 fis depending on the power, and a dark 
NEP ~ 2 x 10 -19 W/\/Hz at the lowest probe power. 
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Quantifying P a b s from the reported results is difficult 
without knowing details of the embedding circuit. How- 
ever our calculations shown in Figs.EJa), (b) and (d) in- 
dicate these densities, temperatures, and lifetimes would 
arise for P a b s in the range 0.1 — 1 aW//im 3 . For the 
same absorbed powers, Fig. |8] indicates a limiting dark 
NEP of 2 - 3 x 10- 19 W/VHz. The agreement with our 
calculations, without free parameters, is extremely satis- 
factory whilst suggesting that the approach that we have 



described has merit. 

In future work we will incorporate a pair-breaking 
source in the kinetic equations in addition to the probe 
signal. We will also investigate the effect of the probe 
frequency on the driven solutions, its effect on the achiev- 
able NEP, and the scaling of 772A with material param- 
eters. It will also be possible to model the detection of 
sub-gap photons < 2A using a driven resonator. 
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